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Coarse grained molecular dynamics simulations are presented in which the sensitivity of the ice nucleation 
rate to the hydrophilicity of a graphene nanoflake is investigated. We find that an optimal interaction strength 
for promoting ice nucleation exists, which coincides with that found previously for an FCC (111) surface. We 
further investigate the role that the layering of interfacial water plays in heterogeneous ice nucleation, and 
demonstrate that the extent of layering is not a good indicator of ice nucleating ability for all surfaces. Our 
results suggest that to be an efficient ice nucleating agent, a surface should not bind water too strongly if it 
is able to accommodate high coverages of water. 


I. INTRODUCTION 

As liquid water is cooled below its melting point it 
crystallizes to solid ice. This familiar yet important pro¬ 
cess is not fully understood, especially at the molecular 
level. It is known that pure water can exist in the liquid 
state far below 0°C and that the reason we see ice for¬ 
mation at temperatures above approximately — 35° C is 
due to the presence of impurity particles! 1 This is known 
as heterogeneous nucleation. It is also known that dif¬ 
ferent impurity particles aid ice formation with different 
efficiencies, for example, feldspar mineral particles have 
been found to be better ice nucleating agents than clay 
mineral particlesP What is severely lacking, however, is 
a comprehensive understanding of heterogeneous ice nu¬ 
cleation: we simply do not understand which properties 
of a material affect its ability to nucleate ice. Given the 
ubiquity of ice formation, and its important role in the 
atmospheric, geological and biological sciences, as well 
as the problems it can cause in the food, transport and 
energy industries, acquiring a full understanding of het¬ 
erogeneous ice nucleation remains a major challenge in 
urgent need of address P 

Experimentally, ice nucleation remains a challenge to 
study as it occurs on small time- and length-scales, 
and computer simulation therefore provides a useful tool 
when investigating both homogeneous 4 ^ and hetero¬ 
geneous ice nucleation. 10 ^ Recent studies have used 
molecular dynamics simulation in combination with 
coarse grained models to probe the mechanisms of het¬ 
erogeneous ice nucleation. In particular, Lupi et al. 13 14 
have investigated the effect of graphitic surfaces on ice 
nucleation, and have found that the extent of layering 
of interfacial water correlates with the freezing tempera¬ 
ture; in our first paper in the series,^ on the other hand, 
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we focused on how the hydrophilicity of an hexagonal sur¬ 
face that acts as a template for the basal face of ice affects 
the rate, and found that an optimal interaction strength 
between the water and the surface exists. In this second 
article, we present further results from simulations of ice 
nucleation in the presence of ‘graphitic’ surfaces. Unlike 
Ref. [HI where the primary aim was to understand ice 
nucleation on soot particles, here we are not attempting 
to model actual graphitic surfaces. Rather, in this study 
we wish to exploit the smoothness of the potential ex¬ 
perienced by the water molecules at such surfaces, and 
compare to the results obtained in the first paper in this 
serieswhere the hexagonal surface under investigation 
presented distinct adsorption sites for the interfacial wa¬ 
ter molecules. We wish to emphasize that we are using 
simplified model surfaces in order to understand possi¬ 
ble general trends that may underlie heterogeneous ice 
nucleation and we therefore probe a far greater range of 
hydrophilicities of these ‘graphitic’ surfaces than previ¬ 
ously considered by Lupi et al. 

The aim of the first paper in this series was to demon¬ 
strate that, by understanding the molecular mechanism 
by which a surface facilitates ice formation, we could 
manipulate the surface to exert a degree of control over 
the rate. The primary purpose of this second article is 
to discuss the results of the first paper in the broader 
context of previous simulations on heterogeneous ice nu¬ 
cleation (in partic ular with respect to the recent work 
of Lupi et a/. 13 * 14 !). In what follows, we will find that 
the graphitic surfaces also exhibit an optimal interaction 
strength with water for promoting ice nucleation, which 
coincides with the optimal interaction strength found for 
the hexagonal surfaces presented in the first paper in this 
seriesWe will al so see that the previously suggested 
layering mechanism 13 14 requires slight modification to 
be applied to strongly adsorbing surfaces and that the 
in-plane structure of the interfacial water molecules can 
affect the layering mechanism. This suggests that the 
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layering mechanism cannot be used to explain the ice nu¬ 
cleating ability of surfaces in general. Finally, we discuss 
the origin of the observed optimal interaction strength 
and suggest a rule-of-thumb for relating the surface hy- 
drophilicity to the ice nucleating efficiency. 


II. METHODS 

A. Systems and force fields 

We have investigated heterogeneous ice nucleation in 
the presence of a rigid graphene nanoflake (GNF) of vary¬ 
ing hydrophilicity, which is totally immersed in water 
as shown in Fig. [l] In this context, an increase in hy- 
drophilicity is synonymous with an increase in the in¬ 
teraction strength between a water molecule and the 
surface. The interaction of water with the GNF was 
modeled using the two-body part of the Stillinger-Weber 
(SW) potential in the same manner as Lupi et a/P 3 ^ 
The GNF consisted of 217 carbon atoms, with a carbon- 
carbon bond-length of 0.142 nm (perfect edges were as¬ 
sumed). The diameter of the GNF was approximately 
2.5 nm to enable direct comparison to the results pre¬ 
sented in the first paper in this series. As in Refs. Il3l 
and m we used asw = 0.32 nm to define the range of 
the water-carbon interaction (the functional form of the 
SW/mW potential is given elsewhere—). The interac¬ 
tion strength was tuned by varying e$w- We note here 
that for the graphitic surfaces in Ref. [I4j values in the 
range 0.12 < esw < 0.2kcal/mol were investigated; as 
this work is concerned with trying to obtain general un¬ 
derstanding rather than modeling a specific system, we 
have broadened this range to 0.06 < esw < 1.5 kcal/mol. 
The total energy after geometry optimization of a single 
water molecule at the center of the GNF was used to de¬ 
fine the water adsorption energy to the surface E a( ^ s (the 
water molecule optimized to a height 0.276 nm above the 
carbon atoms, in the center of a graphene ring). No in¬ 
teraction was defined between the carbon atoms as their 
equations of motion were not integrated. In Ref. [14j it 
was discussed how one can also vary the hydrophilic¬ 
ity of such graphitic surfaces by introducing hydroxyl¬ 
like groups, which leads to different conclusions regard¬ 
ing how the hydrophilicity affects the rate of ice nucle¬ 
ation. The results presented in this article may help un¬ 
derstand this discrepancy, a point that we will return to 
in Sec. |IIIB[ To aid comparison with the work of Lupi et 
aL jiaiiii Table @ shows how F’ads depends upon esw- 

B. Simulation settings 

All simulations were performed using the LAMMPS sim¬ 
ulation package 17 and the coarse grained mW model for 
water. 16 The approximate diameter of an mW water 
molecule is 0.28 nm, as estimated from the radial dis¬ 
tribution function. 16 The velocity Verlet algorithm was 


TABLE I. Dependence of the adsorption energy F a d s on the 
water-carbon interaction strength esw- 


esw (kcal/mol) 

F a ds (kcal/mol) 

0.06 

0.800 

0.13 

1.734 

0.21 

2.801 

0.29 

2.868 

0.37 

4.935 

0.45 

6.002 

0.56 

7.469 

0.67 

8.936 

0.72 

9.603 

0.80 

10.669 

0.88 

11.736 

1.00 

13.337 

1.12 

14.937 

1.31 

17.471 

1.50 

20.005 


used to propagate the equations of motion of the water 
molecules, using a 10 fs time step. In all simulations, 2944 
mW molecules were used and periodic boundary condi¬ 
tions were applied in all three dimensions. As discussed 
in the first paper in this series,^ 2 this system size is suffi¬ 
ciently large that finite size effects do not pose a serious 
problem at this temperature. This includes both increas¬ 
ing the number of water molecules and using a slab geom- 



FIG. 1. (color online) A typical simulation of heterogeneous 
ice nucleation. The GNF is shown by large silver spheres, 
and ice-like molecules are shown by blue lines. The remaining 
liquid-like water molecules are shown by small gray spheres. 
The box shows the boundary of the simulation cell and peri¬ 
odic boundary conditions are used throughout. 
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etry (our results for the GNF are consisten t wit h those 
obtained with a graphitic slab by Lupi et al 13 14 ). Tem¬ 
perature and pressure were maintained using the Nose- 
Hoover thermostat and barostat (with a chain length of 
10) with relaxation times of 1 ps and 2 ps respectively. A 
100 ns trajectory was first performed at 290 K and 1 bar, 
from which initial configurations were drawn (different 
initial configurations were separated by at least 5 ns in 
the high temperature trajectory). At the start of the 
nucleation simulations, velocities for the water molecules 
were drawn randomly from a Maxwell-Boltzmann dis¬ 
tribution to give an initial temperature of 205 K. Sim¬ 
ulations were stopped after 500 ns if nucleation did not 
occur. To detect ‘ice-like’ molecules, we have used the 
CHILL algorithm of Moore et al Rates were extracted 
in the same manner as in the first paper in the seried^ 
and are directly comparable, since we have used the same 
simulation protocol. 


C. Analysis 


For each value of £? a ds> an extra simulation of 10 ns was 
performed at 215 K and lbar (following a 1 ns equilibra¬ 
tion period from a 290K configuration). Similarly, a set 
of 10 ns simulations were performed at 225 K and 1 bar 
for the face centered cubic nanoparticle investigated in 
the first paper in this series.^ (We refer to this nanopar¬ 
ticle as the ‘FCC-111 NP’). These higher temperature 
simulations were performed such that sufficient statistics 
in the liquid state could be obtained over the full range 
of F’ads (i.e. to avoid crystallization over a 10 ns inter¬ 
val). We wish to emphasize that all simulations used to 
calculate the nucleation rate were obtained at 205 K. 

The layering of inter facial water was computed as: 


L = 


L 


^bulk 


dz 


P( Z ) 

Pbulk 


(i) 


where p(z) is the local water density at a height z above 
the surface (see Fig. [2]), Pbuik is the density of bulk liq¬ 
uid water at 215K (or 225K) and lbar (also obtained 
from a 10 ns simulation) and Zbuik = 1.8 nm is a height 
at which p(z) —)► Pbuik- We note that, where comparison 
could be made with the simulations at 205 K, the value 
of L appears to be rather insensitive to the temperature 
(differences are less than 1 unit) - the effect of changing 
-£/ a ds is by far the more dominant effect. The integration 
was performed using the Trapezium rule (Simpson’s rule 
was also used, with a maximum discrepancy between the 
two methods of 3%, and agreement generally within 1%). 
In computing p(z), only water molecules in the column 
above the surface were considered (the radius of the col¬ 
umn was 1.25 nm above both the GNF and the FCC-111 
NP). 

The probability density P(x,y) of water molecules in 
the plane of the surface was computed for the water 
molecules in the contact and second layers above the 



FIG. 2. (color online) Density profile p(z) of water above the 
surface of the (a) GNF at 215 K and (b) the FCC-111 NP at 
225 K, for different values of EMs/AiAap (where AF vap is the 
heat of vaporization of liquid mW water at 298 K). At both 
surfaces water forms layers, with the intensity and sharpness 
of the layers increasing with E a d s - The dark gray shaded 
region indicates water molecules defined as belonging to the 
first layer, and the light gray shaded region wate r mo lecules 
defined as belonging to the second layer (see Sec. IIC). 


surface. A water molecule was defined as being in the 
contact layer if 0 < z < 0.45 nm and in the second layer 
if 0.45 < z < 0.8 nm as shown in Fig. [2](a). A sim¬ 
ilar analysis was also performed for the FCC-111 NP, 
with water molecules defined as being in the contact 
layer if 0 < z < 0.35 nm and in the second layer if 
0.35 < z < 0.7 nm, as shown in Fig. [2](b). 
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III. RESULTS AND DISCUSSION 

A. Nucleation rates and the role of interfacial layering 

In Fig. [3] we show how the nucleation rate R varies 
with the hydrophilicity of the GNF. Specifically, we have 
plotted log-Lo(_R/-^hom) against R a ds/AR va p, where Rhom 
is the homogeneous nucleation rate and AF vap is the 
heat of vaporization of bulk mW water (10.65 kcal/mol 
at 298K)P= 6 We can clearly see that for the weakest in¬ 
teraction strength, the GNF has little effect on the nucle¬ 
ation rate. As R a d s increases, the rate rapidly increases 
to reach a maximum at R a d s /AR vap & 0.3—0.4 that is 
approximately a factor 25 faster than homogeneous nu¬ 
cleation. We note here that a factor 25 increase in the 
rate would appear small when compared to experimental 
values, which often span many orders of magnitude. This 
is due to the fact that we are operating at low temper¬ 
atures so that we can directly compare to homogeneous 
nucleation. The effects of heterogeneous ice nucleation 
will become more pronounced at higher temperatures. 
Upon increasing Rads further, the rate steadily drops un¬ 
til Rads/ARvap ~ 1.0 when the rate begins to steadily 
increase again. For the most strongly interacting GNF 
investigated, the rate is a little over 10 times that of 
homogeneous nucleation. We also show the results from 
Ref. [12] of the nucleation rate in the presence of the FCC- 
111 NP. Both the GNF and the FCC-111 NP exhibit a 
maximum in rate at R a d s /AR vap ~ 0.3—0.4, but differ 
in that at R a d s /AR vap ~ 1-0, the FCC-111 NP crosses 
from promoting to inhibiting rather than exhibiting a 
local minimum like the GNF. This crossover from pro¬ 
moting to inhibiting can be explained by an excess of 
favorable adsorption sites at the FCC-111 NP (shown in 
Fig. El and discussed in detail in Ref. EE2). 

To explain the ice nucleating ability of the graphitic 
surfaces such as those considered here, Lupi et al. found 
that the layering of interfacial water L (see EclTI) cor¬ 
relates with the ice nucleating abilityP3 In Fig. [5 (a), we 
show the dependence of L on R a d s for the GNF. Un¬ 
surprisingly, L increases monotonically with R a ds and we 
therefore cannot explain the observed non-monotonic de¬ 
pendence of R on R a ds seen in Fig. [3] simply by the ex¬ 
tent of layering. Figs. [6](a) and [6](b) provide some in¬ 
sight into why this is the case, where we show typical 
structures of water in contact with the GNF at 215 K for 
R a ds/ARv a p « 0.25 and R ads /AR vap « 1.9, respectively. 
For values of R a d s that yield the highest rates, such as in 
Fig ,[6](a), water forms structures in the contact layer that 
resemble the hexagonal structure of ice. Indeed, when ice 
formation is observed at 205 K, it appears to be driven by 
the formation of such hexagonal patch es in the contact 
layer, consistent with previous studies.^ 3 331 In the case 
of the strongly adsorbing GNF shown in Fig. [6](b), it is 
clear that the number of water molecules in contact with 
the GNF has increased and that, rather than an hexag¬ 
onal structure similar to ice, a structure consisting pre¬ 
dominantly of smaller membered rings is now observed. 
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FIG. 3. (color online) Dependence of the heterogeneous nu¬ 
cleation rate on surface hydrophilicity. As R a ds increases so 
too does the hydrophilicity. The homogeneous and FCC-111 
data are taken from Ref. FIR! Like the FCC-111 NP, the GNF 
also exhibits a maximum rate at R a ds/AR V ap ~ 0.3—0.4, 
but in contrast, exhibits a local minimum in the rate at 
Rads/AR va p » 1.0. 


This structure bears a strong resemblance to those seen 
in confined water layers under high pressures 19 20 and can 
be understood as water maximizing its interaction to the 
surface with only a slight cost in hydrogen bond energyP^I 
When ice formation is observed at this surface, it does 
so in the layers of water above the contact layer, with 
the structure in the contact layer remaining unchanged. 



FIG. 4. (color online) Typical structures that form in the 
contact layer at the FCC-111 NP at 225 K. The FCC-111 NP 
is shown in silver and the water molecules in blue, (a) FCC- 
111 NP with R a ds/ARvap ~ 0.3. An hexagonal overlayer, 
which is commensurate with the surface, that resembles ice 
is observed and facilitates ice formation at 205 K. Unoccu¬ 
pied, or ‘excess’, adsorption sites (highlighted by yellow cir¬ 
cles) are present when this structure forms, (b) FCC-111 NP 
with R a ds/AR V ap ~ 1.9. For this stronger interaction with 
the surface, the excess sites are occupied and the hexagonal 
structure resembling ice is no longer observed. 
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FIG. 5. (color online) Dependence of the extent of layering 
of interfacial water on E a d s - The black squares show the lay¬ 
ering over the whole density profile L as defined by Eq. |TJ 
whereas the red circles show the layering excluding the con¬ 
tact layer L* as defined by Eq. [ 5 ] (a) Result for the GNF at 
215 K. Both L and L * increase monotonically with E a d s , but 
L does so much more rapidly. The inset contains the same 
data but with a smaller scale for the y- axis, which shows 
more clearly that L* > 3 only once F a ds/Ai7 V ap ^ 0.7 (when 
E a ds/Ai7vap ~ 0.1, L 3 and nucleation is not enhanced), 
(b) Results from the FCC-111 NP at 225 K. Although L is 
much greater at the FCC-111 NP than at the GNF, the val¬ 
ues of L* are comparable. 


Thus for high values of E a d s , the contact layer is ‘inactive’ 
with respect to nucleation. 

The observation that the contact layer becomes inac¬ 
tive to ice nucleation for strong adsorption energies is 
enough to understand why we begin to see a decrease in 
the nucleation rate beyond E a d s /AH vap ~ 0.3—0.4. To 
explain the increase in rate beyond ^ a ds/Ai7v ap ~ 1.0, 



FIG. 6. (color online) Typical structures of water in the 
contact layer at the GNF at 215 K. The GNF is shown 
in silver and the water molecules in blue, (a) GNF with 
F a ds/Ai7vap ~ 0.25. Patches of ice-like hexagons readily form 
in the contact layer and facilitate ice formation at 205 K. (b) 
GNF with E a ds/AH vap 1.9. At this more strongly adsorb¬ 
ing surface, the structure in the contact layer consists predom¬ 
inantly of smaller membered rings. This persists even after 
ice nucleation at 205 K, which at this interaction strength, 
occurs in the water layers above. 


however, requires further analysis. To this end, we have 
computed the layering of interfacial water with contribu¬ 
tions from the first layer excluded: 

rz bulk 

L* = / dz 

j Z 0 

where, at the GNF, zq = 0.45 nm. In Fig. |5](a) we can 
see that L*, like L, also increases monotonically with 
-E/ a ds •> but much more slowly. We can also see that the 
value of L* at E a d s /AEL^^ ~ 1.9 is similar to the value 
of L at E a ds/AH vap ~ 0.2 and from Fig. [31 that these 
two adsorption energies yield similar rates (both approx¬ 
imately a factor 10 faster than homogeneous nucleation). 
It therefore seems that beyond E a d s /AH vap « 1.0, the 
extent of layering in the second layer of water and above 
becomes sufficient to promote ice nucleation. This can 
be seen more clearly in the inset of Fig. [5](a); bearing in 
mind that the most weakly interacting GNF yields L « 3 
and does not promote ice nucleation, we can see that L* 
only begins to exceed this value for E a d s / AH vap >0.7. 


p{z) 


- 1 


Pbulk 


( 2 ) 


B. The layering mechanism depends upon the in-plane 
structure of water 

Section |III A| provides strong evidence in support of 
the layering mechanism, albeit with a sl ight modifica¬ 
tion to what was originally proposed. b^lMI Conceptually, 
the layering mechanism is appealing and can perhaps be 
understood in terms of reducing the entropic barrier to 
nucleation: if water molecules are restricted to motion 
in a particular plane (with a density acceptable for ice 
nucleation), then the space that they can explore is effec¬ 
tively reduced by one dimension relative to the bulk liq- 
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uid, which subsequently reduces the number of possible 
configurations that the water molecules can explore. This 
argument, however, implicitly assumes that the effective 
potential experienced by the water molecules within a 
layer is uniform. 

For the graphitic surfaces investigated in this study, 
the assumption of a uniform effective potential within the 
layers is likely to be reasonable. Now consider the FCC- 
111 NP with F’ads/Ai^vap ^ 1.9 which, as can be seen in 
Fig-! does not promote ice nucleation. The GNF with 
similar E a d s , on the other hand, enhances ice nucleation 
by a factor ~10 relative to homogeneous nucleation. Fur¬ 
thermore, the layering (excluding contributions from the 
contact layer) above both of these surfaces is similar, with 


L* « 6 for the FCC-111 NP and L* « 5 for the GNF, 
as shown in Fig. [5](b) (a value of zq = 0.35 nm is used in 
Eq. [2] for the FCC-111 NP). The difference in rates can 
be understood in terms of in-plane structure, as seen in 
Fig. [7, where we show — ln[P(x, y)\, where P(x,y) is the 
probability density of water molecules in the plane of the 
surface at 215 K, both for the contac t laye r and the sec¬ 
ond layer above the surface (see Sec. IIC). At the FCC- 
111 NP, the water molecules bind at distinct adsorption 
sites (see Fig. 0(b)) and do not diffuse over the 10 ns 
timescale of the simulation. Importantly, this impacts 
upon the structure of the water molecules in the second 
layer (see Fig. 0(d)), which tend to be found directly 
above those in the contact layer. In contrast, the water 
molecules in contact with the GNF (see Fig. 0(a)) do 
not adsorb at particular adsorption sites and are not im¬ 
mobile like at the FCC-111 NP, resulting in a smearing- 
out of P(x,y). Accordingly, P(x,y) for the second layer 
above the GNF (see Fig.[7](c)) is much smoother than at 
the FCC-111 NP. The smoothness of P{x,y) for the sec¬ 
ond layer above the GNF means that the water molecules 
can rearrange to form ice-like structures, whereas the cor¬ 
rugation in P(x, y) for the second layer above the FCC- 
111 NP appears to frustrate the water molecules, hinder¬ 
ing ice formation. 


In Refs. fl3l and [T4l the general applicability of the lay¬ 
ering mechanism was left as an open question. By vastly 
broadening the range of hydrophilicities investigated and 
monitoring the response of the ice nucleation rate in the 
presence of two model surfaces, we are able to elucidate 
when the layering mechanism may be important. The 
results of our simulations show that the extent of layer¬ 
ing can be important for ice nucleation, but that if the 
coverage of water at the surface is too high, then the con¬ 
tributions of the contact layer to the layering should be 
omitted. Furthermore, the importance of layering is de¬ 
pendent upon the water molecules experiencing a rather 
uniform effective potential within the layers. The extent 
of layering is therefore not a universal descriptor for the 
ice nucleating ability of surfaces. 

We finish this section with a comment regarding the 
manner by which the hydrophilicity of the GNF has 
been tuned. In this study, this has been achieved by 
uniformly changing the value of esw- In Ref. na how¬ 



FIG. 7. (color online) In-plane distribution of water mole cules 
above the surface, plotted as — ln[P(;c, y)\ (see Sec. [net with 
Pads/Aifvap ~ 1.9 at 215 K. (a) and (b) show the contact 
layer at the GNF and FCC-111 NP, respectively, (c) and (d) 
show the second layer above the GNF and FCC-111 NP, re¬ 
spectively. Unlike at the GNF, the water molecules in the 
contact layer with the FCC-111 NP are bound at specific ad¬ 
sorption sites and do not diffuse over the timescale of the 
simulation (10 ns). The water molecules in the second layer 
above the FCC-111 NP consequently exhibit greater structure 
than those in the second layer above the GNF. 


ever, Lupi et al. also investigated the effect of introduc¬ 
ing hydroxyl-like groups (with higher concentrations of 
hydroxyl-like groups corresponding to more hydrophilic 
surfaces). They found that increasing the hydrophilic¬ 
ity in such a manner was in fact detrimental to the 
rate, while changing esw over the range 0.12-0.2 keal/mol 
enhanced ice nucleation, consistent with our findings 
(see Table I|. The results presented in this section may 
reconcile the discrepancy between the two approaches: 
by uniformly increasing esw over this relatively narrow 
range, ice nucleation is enhanced due to an increase in 
L as the water molecules are still moving in a relatively 
smooth effective potential; the introduction of hydroxyl¬ 
like groups, on the other hand, is likely to localize water 
molecules at certain positions at the surface, which may 
cause a similar frustration to that described at the FCC- 
111 NP, if the spatial arrangement of hydoxy 1-like groups 
is not conducive to ice nucleation. 


C. The effect of surface hydrophilicity on heterogeneous 
ice nucleation 

In this section, we will discuss the observation that 
the GNF and the FCC-111 NP both have optimal in¬ 
teraction strengths at E a( \ s /AH vap « 0.3—0.4 in more 
detail. If the interaction between the water and the sur¬ 
face is too weak, the induced structural differences from 
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the bulk liquid are not significant enough to promote ice 
nucleation at either the GNF or the FCC-111 NP. What 
is more intriguing is why the ice nucleation rate at both 
surfaces decreases beyond E a a s /AH vap « 0.3—0.4. De¬ 
spite the differences in surface topography (the GNF can 
be considered a smooth surface whereas the FCC-111 NP 
presents distinct adsorption sites), both surfaces share 
a common feature: they can both accommodate water 
coverages that are higher than that when ice forms at 
the surface. This has been demonstrated qualitatively 
in Figs. [4] and i Fig-i provides quantitative evidence 
for this statement, where we show the lateral density of 
water molecules a in the contact and second layers: 



whe re th e integral runs over the layer of interest (see 
Sec. IIC). Note that a can also be computed explicitly 
by counting the average number of water molecules in 
a given layer: such an approach also gives information 
regarding the fluctuations of a. For the GNF, shown 
in Fig. [8](a), we can see that a for the contact layer 
steadily increases with E a ^ s whereas for the second layer, 
although increasing slightly initially, a is essentially con¬ 
stant. In terms of layering, this means L* is increasing 
primarily through a narrowing of the second peak in p(z), 
whereas L also has a contribution from an increased num¬ 
ber of water molecules at the surface. To a lesser extent 
the same is true for the FCC-111 NP, shown in Fig.[8](b), 
although some variation in a for the second layer is also 
observed. 


We have also marked in Figs. [8](a) and (b) the water 
coverage of ice that forms at the surface for the GNF and 
the FCC-111 NP, respectively.^ We can clearly see that 
for F^ads/AiT’vap « 0.3—0.4, the value of a in the contact 
layer is comparable to that of ice that forms at the sur¬ 
face. Thus, as F'ads increases further, it becomes increas¬ 
ingly favorable for water to adsorb to the surface, to the 
detriment of ice nucleation occurring in the contact layer. 
We therefore suggest a rule-of-thumb for the role of sur¬ 
face hydrophilicity in ice nucleation: for surfaces that can 
accommodate water coverages higher than that required 
by ice, binding to the surface should not be too strong, 
with optimal adsorption energies approximately 30-^0% 
the heat of vaporization of liquid water. We must stress 
that the importance of surface hydrophilicity will depend 
upon the water coverage that the surface can accommo¬ 
date. As we have shown in the first paper in this series,^ 
surfaces with E a ^ s /AH vap ^ 0.4 can exhibit excellent ice 
nucleating efficiency, provided that the coverage of water 
at the surface does not exceed that of ice. For example, 
surfaces that resemble the surface of ice itself may also 
favor more open structures, such that this rule-of-thumb 
cannot be applied. 22 
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FIG. 8. (color online) Dependence of the lateral density of 
water molecules a on F a d s for (a) the GNF at 215 K and (b) 
the FCC-111 NP at 225 K. The blue dashed lines indicate 
the water coverage of ice that forms at both surfaces (the 
blue shaded area in (a) indicates the standard deviation of a 
sample average, taken over the range 0.16 < F a d s /Ai7 V ap < 
0.56). The black and red shaded areas indicate the standard 
deviation in the measured values of a for the contact and 
second layers, respectively. When F a d s /AI4vap ~ 0.3—0.4, the 
coverage in the contact layer is close to the water coverage of 
ice that forms. 


IV. CONCLUSIONS 

We have presented computer simulations of hetero¬ 
geneous ice nucleation in the presence of a graphene 
nanoflake immersed in water and investigated how its 
hydrophilicity affects the nucleation rate. The results of 
our simulations in part support thepreviously proposed 
layering mechanism of Lupi et al. j 13 l 14 although we have 
seen that for strongly adsorbing surfaces, the increased 
layering, due to a higher coverage of water molecules, is 










detrimental to ice nucleation. Under such conditions, by 
excluding the contribution of the water molecules in con¬ 
tact with the surface, the extent of layering is, however, 
still found to correlate with the heterogeneous nucleation 
rate. It has also been demonstrated that the layering 
mechanism is not universal, as surfaces that exhibit simi¬ 
lar degrees of interfacial layering can yield vastly different 
rates. We attribute this finding to the extent that the sur¬ 
face affects the in-plane structure of the water molecules: 
for surfaces where the water molecules move in a smooth 
effective potential, like the graphitic surfaces investigated 
in this article, the extent of layering describes the nucle¬ 
ation rate well; for surfaces that present distinct adsorp¬ 
tion sites, such as the FCC (111) surface investigated in 
Ref. [TH the induced structure can frustrate ice nucleation 
not only in the contact layer, but also in the layer above. 

We have observed that an optimal interaction strength 
between the graphene nanoflake and water for ice nu¬ 
cleation exists, and that this coincides with the optimal 
interaction strength found for the FCC (111) surface in¬ 
vestigated previously.^ This behavior has been rational¬ 
ized by noting that both surfaces are able to accommo¬ 
date coverages of water that are higher than that when 
ice forms at these surfaces. We have proposed a rule- 
of-thumb, which states that in order to nucleate ice ef¬ 
ficiently, a surface should not bind water too strongly 
if it can accommodate high coverages of water. Such in¬ 
sight may prove useful when trying to predict a material’s 
ice nucleating ability, especially as the coverage of liquid 
water should be obtainable through e.g. surface X-ray 
diffraction experiments. 
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